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1 Theory 

1.1 Introduction 

This discussion serves as an introduction to the use of Monte Carlo simulations 
as a useful way to evaluate the observables of a ferromagnet. Key background 
is given about the relevance and effectiveness of this stochastic approach and in 
particular the applicability of the Metropolis-Hastings algorithm. Importantly 
the potentially devastating effects of spontaneous magnetization are highlighted 
and a means to avert this is examined. 

An Ising model is introduced and used to investigate the properties of a 
two dimensional ferromagnet with respect to its magnetization and energy at 
varying temperatures. The observables are calculated and a phase transition 
at a critical temperature is also illustrated and evaluated. Lastly a finite size 
scaling analysis is underatken to determine the critical exponents and the Curie 
temperature is calculated using a ratio of cumulants with differing lattice sizes. 
The results obtained from the simulation are compared to exact calculations to 
endorse the validity of this numerical process. A copy of the code used, written 
in CH — h, is enclosed and is freely available for use and modification under the 
General Public License. 

1.2 Background 

In most ordinary materials the associated magnetic dipoles of the atoms have a 
random orientation. In effect this non-specific distribution results in no overall 
macroscopic magnetic moment. However in certain cases, such as iron, a mag- 
netic moment is produced as a result of a preferred alignment of the atomic 
spins. 

This phenomenon is based on two fundamental principles, namely energy 
minimization and entropy maximization. These are competing principles and 
are important in moderating the overall effect. Temperature is the mediator 
between these opposing elements and ultimately determines which will be more 
dominant. 

The relative importance of the energy minimization and entropy maximiza- 
tion is governed in nature by a specific probability 

P(a)=exp^^. (1) 
which is illustrated in Figure [Hand is known as the Gibbs distribution. 

1.3 Model 

A key ingredient in the understanding of this theory is the spin and associated 
magnetic moment. Knowing that spin is a quantum mechanical phenomenon it 
is easy to anticipate that a complete and thorough exposition of the problem 
would require quantum rules of spin and angular momentum. These factors 
prove to be unnecessary complications. 

We thus introduce a model to attain useful results. The idea central to a 
model is to simplify the complexity of the problem to such a degree that it is 
mathematically tractable to deal with while retaining the essential physics of 
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Boltzmann Probability Distrubution vs Energy vs Temprature 




Figure 1: This figure shows the Boltzmann probabiHty distribution as a land- 
scape for varying Energy (E) and Temperature (T). 



the system. The Ising Model does this very effectively and even allows for a 
good conceptual understanding. 

t t t I 
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Figure 2: Two Dimensional lattice illustration of an Ising Model. The up and 
down arrows represent a postive and a negative spin respectively. 

The Ising Model considers the problem in two dimensions^ and places dipole 
spins at regular lattice points while restricting their spin axis to be either up 
(+y) or down (-y) . The lattice configuration is square with dimensions L and the 
total number of spins equal to N = L x L. In its simplest form the interaction 
range amongst the dipoles is restricted to immediately adjacent sites (nearest 
neighbours). This produces a Hamiltonian for a specific spin site, i, of the form 

Hi = -jJ2siSj (2) 

where the sum j„„ runs over the nearest neighbours of i. The coupling constant 
between nearest neighbours is represented by J while the Sj and sj are the re- 
spective nearest neighbour spins. The nature of the interaction in the model is 
all contained in the sign of the interaction coupling constant J. If J is positive 



^It may also be considered in three dimensions but it should be noted that a one dimensional 
representation does not display phase transition. 
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it would mean that the material has a ferromagnetic nature (parallel align- 
ment) while a negative sign would imply that the material is antiferromagentic 
(favours anti-parallel alignment). J will be taken to be -1-1 in our discussion 
and the values for spins will be -1-1 for spin up and -1 for spin down. A further 
simplification is made in that J/kb is taken to be unity. The relative positioning 
of nearest neighbours of spins is shown in Figure [3] with the darker dot being 
interacted on by its surrounding neighbours. 

ABOVE 
(x;y+l) 

O 
A 



(x-i;y) A J (x+i;y) 



V 

o 

BELOW 
(x; y-1) 

Figure 3: Nearest neighbour coupling. The dark dot, at position (x,y), is being 
interacted upon by its nearest neighbours which are one lattice spacing away 
from it. 



To maximize the interaction of the spins at the edges of the lattice they are 
made to interact with the spins at the geometric opposite edges of the lattice. 
This is referred to as periodic boundary condition (pbc) and can be visualized 
better if we consider the 2d lattice being folded into a 3d torus with spins being 
on the surface of this topological structure. 



Figure 4: An illustration of a three dimensional torus which is repreentative of 
a two dimensional lattice with periodic boundary conditions. 



1.4 Computational Problems 

With the help of an Ising Model we can proceed with the anticipation of achiev- 
ing solutions for observables. If the energy of each possible state of the system 
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is specified, then the Boltzmann distribution function, equation ([T]), gives the 
probabihty for the system to be in each possible state (at a given temperature) 
and therefore macroscopic quantities of interest can be calculated by doing prob- 
ability summing. This can be illustrated by using magnetization and energy as 
examples. For any fixed state, a, the magnetization is proportional to the 'ex- 
cess' number of spins pointing up or down while the energy is given by the 
Hamiltonian 

M{a) = Nup{a) - Ndown{a) (3) 
The expected value for M is given by 

(A/)-^M(a)P(a), (4) 

a 

and the expected value for E is given by 

(i?)=5]£;(a)P(a). (5) 

a 

These calculation pose a drastic problem from a practical perspective. Consid- 
ering we have two spin orientations (up & down) and there are N spins which 
implies that there are 2^ different states. As N becomes large it is evident that 
computations become a daunting task if calculated in this manner. 

It may seem a natural suggestion to use a computer simulation to do these 
calculations but by examining equations (jlj and ([5]) more closely it becomes 
apparent that using this procedure would waste as much computing effort on 
calculating an improbable result as it does on a very probable result. Thus a 
better numerical alternative would be to use a simulation to generate data over 
the 'representative states'. These representative states constitute the appropri- 
ate proportions of different states ^. This is a form of biased sampling which 
essentially boils down to satisfying the following condition 

GENERATED FREQUENCYeeACTUAL PROBABILITY, 
(computer) (theory) 

We now examine, in a more formal setting, how to accomplish this objective. 



1.5 Sampling and Averaging 

The thermal average for an observable A{x) is defined in the canonical ensemble 

{A{x))t = I y e-P"<^-^A{x)dx (6) 

where a; is a vector in the phase space and /? — 1/khT. The Partition function, 
Z, is given by 

Z = Je-^"^-Ux 
while the normalized Boltzmann factor is 



^The frequency with which some class of events is encountered in the representative sum 
must be the same as the actual probability for that class. 
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P{x) = ^e-^"(^\ (7) 

This probability gives the actual statistical weight with which the configuration 
X occurs in the thermal equilibrium. We now want to consider the discrete case 
of the formal definitions above. If we are to consider a finite portion of the 
phase space it would produces an average of the form 



= ^':^M ...:r ■ (8) 



If we were to take M ^ cxo in equation ([5]) it would reduce to equation ^ . The 
problem with taking a simple sample of this nature in the phase space is that 
it would not guarantee that the probability distribution is peaked in the region 
under consideration (not representative). Figure [H illustrates this problem. 



P(E) 




<H> E 

N 



Figure 5: Example of a simple sampling producing a Gaussian distribution, 
centered around zero, while the crucial data is peaked outside the sampling 
region. 



It thus makes sense to attempt a smarter sampling technique to include the 
important areas in the phase space. We want a process that selects points, 
xi, with an associated probability, P{xi) in the phase space. Estimating the 
thermal average now for a chosen set, a;;, reduces equation ([5]) to 



(^(-)) = ^':^M ..,:r^ . ■ (9) 



Y.tUe-"^^'^^A{xi)/P{x i) 

The most sensible choice for P{xi) is P{xi) oc e"^*^^'-''^. This construct produces 
the simple arithmetic average for Q by canceling out the Boltzmann factors, 
thxis 

M 

= mE^(^')- (10) 

1=1 
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If we stop to reflect on what we are trying to achieve at this point, we discover 
that we are attempting to reduce a probabiUty distribution at equihbrium of the 
infinite phase space to a representative distribution with a finite set of points 
from the phase space, xi . The question is how do we generate this distribution. 

Metropohs et al. advanced conventional thinking at the time by introducing 
the idea of using a Markov process of successive states xi where each state 
xi^i is constructed from a previous state xi via a suitable transition probability 
W{xi xm-i). In order to implement this idea successfully a detailed balance 
equation has to be imposed, 

Peq{xi)W{xi ^ Xl') = Peq{xi')W{xi> ^ Xi). (11) 

This equation simply prescribes that at equilibrium there is an equal probability 
for xi — > xi' and xi' — > xi. If we now take the ratio of the transition probabilities 
it becomes evident that a move xi — ^ xi' and the inverse move Xf — > xi is only 
dependent on the energy change SH = H{xi) — H{xi). 

W{xi ^Xv) _ 



W{xv ^ xi) 



e-'"^ (12) 



This doesn't however help to specify W{xi xii) uniquely. In order to do this 
we introduce 

r e~*^^ if < , 
W{xi ^ XI.) = \ (13) 
I 1 otherwise , 

It can be shown that using this transition probability W{xi a;;+i) the distri- 
bution P{xi) of states generated by the Markov process tend to the equilibrium 
distribution as M — > oo. Thus the construct holds and approximates the theory 
with an increasing degree of accuracy as we consider a larger set of points, {xi}, 
in the phase space. 

The changes in the probability distribution over time are governed by the 
Markovian Master equation 

= - ^ W{x, ^x,) + Y: Wix, - X,). (14) 

x' x' 

In the thermal limit where P{xi) — Peq{xi) the detailed balance equation, (fTTj) . 
which was imposed, comes into play resulting in dP{x,t)/dt = as we would 
expect. Furthermore since we are considering a finite system it is logical to 
conclude that the system is ergodic^. 

{A{t)) = 1 A{t)dt (15) 

This relation reduces in a similar fashion to the arithmetic average previously 
discussed if we consider the number of Monte Carlo steps (mcs) as a units 
of 'time'. We thus are confronted by the question of whether the system is 



^Attribute of a behavior that involves only equilibrium states and whose transition prob- 
abilities either are unvarying or follow a definite cycle. In statistics, ergodicity is called 
stationarity and tested by comparing the transition probabilities or different parts of a longer 
sequence of events. 
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ergodic in order for the time averaging to be equivalent to the canonical ensemble 
average. This condition is thus forced upon the system if we consider the mcs 
as a measure of time. 

M 

(^W) = mE^(^W) (16) 

Thus the Metropolis sampling can be interpreted as time averaging along a 
stochastic trajectory in phase space, controlled by the master equation of 
the system. 



1.6 Monte Carlo Method 

The fact that we can make a successful stochastic interpretation of the sampling 
procedure proves to be pivotal in allowing us to introduce the Monte Carlo 
method as a technique to solve for the observables of the Ising Model. 

A Monte Carlo calculation is defined, fundamentally, as explicitly using ran- 
dom variates and a stochastic process is characterized by a random process that 
develops with time. The Monte Carlo method thus lends itself very naturally 
to simulating systems in which stochastic processes occur. From what we have 
established with regards to the sampling being analogous to a time averaging 
along a stochastic trajectory in the phase space it is possible to simulated this 
process by using the Monte Carlo method. The algorithm used is design around 
the principle of the Metropolis sampling. 

From an implementation point of view the issue of the random number lies at 
the heart of this process and its success depends on the fact that the generated 
number is truly random. 'Numerical Recipes in C, [9], deals with this issue 
extensively and the random number generator 'Ranl.c' from this text was used 
in the proceeding simulation. This association with randomness is also the origin 
of the name of the simulation since the glamorous location of Monte Carlo is 
synonymous with luck and chance. 



1.7 Calculation of observables 

The observables of particular interest are (E), (i?^), (A/), (|M|) and (A'P). 
These are calculated in the following way, 

1 ^ 

a 

similarly (|Af|) and {M'^) are calculated using the above equation. 
To calculate energy we use the Hamiltonian given in equation ([2]) , 

N 1 ^ 

(E) = 2^ ) - 2^ -JT.J2'^'^ ) (18) 

i i jnn 

the factor of a half is introduced in order to account for the spins being counted 
twice. Equation psp is used in a similar way to determine (i?^). 

At the Curie temperature we expect a marked fluctuation in these quan- 
tities. A good candidate to illustrate this fluctuation would be the variance 
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(AA)^ = (A^) — (A)'^. This leads us to the logical conclusion of calculating the 
heat capacity, C, and the susceptibility, x- 



_ (AE)^ _ (E^) ~ {Ey 

dM (AM)2 (aP) - (M)2 



X 



dT 



hT 



hT 



(19) 
(20) 



A cumulant is also calculated. This will be used to ultimately determine the 
Curie temperature. 



(21) 



1.8 Metropolis Algorithm 

The algorithm implemented for this simulation is the Metropolis Algorithm [8] . 
The steps executed in the program are best summarized in a flowchart. From 
the flowchart. Figure [HI it is possible to attain a better conceptual feel for what 
the algorithm is attempting to achieve. 



Intialize 
Lattice 





End Program 











/ Ouput 
data / 




Average 
observables 











Figure 6: Metropolis Flowchart 
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• In the first step the lattice is INITIALIZED to a starting configuration. 

This may cither be a homogencoiis or random configuration. A random 
configuration has a benefit in that it uses less computing time to reach a 
equilibrated configuration with the associated heat bath. 

• In the following PROCESS the random number generator is used to se- 
lect a position on the lattice by producing a uniformly generated number 

between 1 and N. 

• A DECISION is then made whether the change in energy of flipping that 
particular spin selected is lower than zero. This is in accordance with the 
principle of energy minimization. 

— If the change in energy is lower than zero then a PROCESS is invoked 

to flip the spin at the selected site and the associated change in the 
observables that want to be monitored are stored. 

— If the change in energy is higher than zero then a DECISION has to 
be used to establish if the spin is going to be flipped, regardless of the 
higher energy consideration. A random number is generated between 
and 1 and then weighed against the Boltzmann Probability factor. 
If the random number is less than the associated probability, e~^^^ , 
then the spin is flipped (This would allow for the spin to be flipped 
as a result of energy absorbed from the heat bath, as in keeping with 
the principle of entropy maximization) else it is left unchanged in its 
original conflguration. 

• The above steps are repeated N times and checked at this point in a 
DECISION to determine if the loop is completed. The steps referred to 
here do not include the initialization which is only required once in the 
beginning of the algorithm. 

• Once the N steps are completed a PROCESS is used to add all the pro- 
gressive changes in the lattice configuration together with the original 
configuration in order to produce a new lattice configuration. 

• All these steps are, in turn, contained within a Monte Carlo loop. A 
DECISION is used to see if these steps are completed. 

• Once the Monte Carlo loop is completed the program is left with, what 
amounts to, the sum of all the generated lattices within the N loops. A 
PROCESS is thus employed to average the accumulated change in observ- 
ables over the number of spins and the number of Monte Carlo steps. 

• Lastly this data can be OUTPUT to a file or plot. 

This run through the algorithm produces a set of observables for a specific 
temperature. Considering that we are interested in seeing a phase transition 
with respect to temperature we need to contain this procedure within a temper- 
ature loop in order to produce these observables for a range of temperatures. 

The program that was implemented for this discussion started at a temper- 
ature of T = 5 and progressively stepped down in temperature to T = 0.5 with 
intervals of 5T = 0.1. The different lattice sizes considered where 2x2, 4x4, 
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8x8 and 16 x 16. A million Monte Carlo steps (mcs) per spin where used in 
order to ensure a large sample of data to average over. The physical variables 
calculated were (E), (E^), (|M|) and (M^). 

A problem that occurs after the initialization is that the configuration will, 
more than likely, not be in equilibrium with the heat bath and it will take a few 
Monte Carlo steps to reach a equilibrated state. The results produced during 
this period are called transient and aren't of interest. We thus have to make 
provision for the program to disregard them. This is achieved by doing a run of 
a thousand Monte Carlo steps preceding the data collection of any statistics for 
a given temperature in order to ensure that it has reached a equilibrated state. 
This realization is only significant for the initial configuration and the problem is 
avoided by the program in the future by using small temperature steps and the 
configuration of the lattice at the previous temperature. A very small number 
of mcs are thus required for the system to stabilize its configuration to the new 
temperature. 



2 Results 

2.1 Energy Results 
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Figure 7: This plot shows the differing results of the Energy for varying lattice 
sizes, L X L. 



In Figure [7] the energy per spin as a function of temperature can be seen. 
The curve of the graph becomes more pronounced as the lattice size increases 
but there isn't a marked difference between the L — 8 and L ~ 16 lattices. The 
steep gradient in the larger lattices points towards a possible phase transition 
but isn't clearly illustrated. The energy per spin for higher temperatures is 
relatively high which is in keeping with our expectation of having a random 
configuration while it stabilizes to a E/N = —2J — —2 at low temperatures. 
This indicates that the spins are all aligned in parallel. 
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Specific Heat Capacity per spin (C/N) vs Temperature (T) 
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Figure 8: This plot shows the differing results of Specific Heat Capacity for 
varying lattice sizes, L x L. 



In Figure [8] the specific heat capacity per spin is shown as a function of tem- 
perature. We concluded previously that a divergence would occur at a phase 
transition and thus should be looking for such a divergence on the graph. It is 
however clear that there is no such divergence but merely a progressive steep- 
ening of the peak as the lattice size increases. The point at which the plot is 
peaked should be noted as a possible point of divergence. The reason for not 
explicitly finding a divergence will be discussed in Section [^T^ 



2.2 Magnetization Results 

Figure [5] of the magnetization results shows very beautifully that the shape 
of the gradient becomes more distinct as the lattice size is increased. Fur- 
thermore, as opposed to Figure H there is a far more apparent difference that 
the larger lattices produce in the curves and this illustrates a more apparent 
continuous phase transition. The behaviour of the magnetization at high and 
low temperature are as the theory prescribes (random to stable parallel aligned 
configuration) . 

At this juncture it is prudent to point out that the susceptibility cannot be 
calculated using the ordinary technique in equation ([20]) given in the discussion 
on the calculation of observables. The reason is focused around a subtle fact 
that has drastic implications. To comprehend the problem at work we have to 
consider one of the constraints of our model, namely the finite nature of our 
lattice. This manifests in the fact that spontaneous magnetization can occur 
for a finite sized lattice. In this instance the effect is of particular interest below 
the critical temperature. 

This can be illustrated by considering the following example of collected 
data in Figure [TOl This data is taken at a temperature that is considerably 
less than the Curie temperature and we would thus expect it to have a stable 
nature and yet it clearly displays a fluctuation that is uncharacteristic, resulting 
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Absolute Magnetisation per spin (<|I\/1|>/N) vs Temperature (T) 
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Figure 9: This plot shows the differing resuhs of the Magnetization for varying 
lattice sizes, L x L. 

in a complete flip of the magnetization. It has already been highlighted that 
because we are dealing with a limited lattice size there is finite probability for 
this kind of behaviour to take place. This probability is directly proportional 
to the number of mcs used and inversely proportional to the lattice size, this is 
compounded by the preiodic boundary conditions used. 

Figure [TT] schematically depicts this fact. The valley shown linking the two 
peaks of the probability will thus be dropped for lower temperatures and bigger 
lattice configurations. It should be noted that even though the probability may 
be less it does always exist and this has to be accounted for in data collection or 
it may corrupt the results. We expect the configuration to be relatively stable 
at the peaks but if its magnetization has slipped down (fluctuated) to the center 
of the valley then it has an equal probability of climbing up either side of the 
peaks, this is the crux of the spontaneous flipping. This aspect of symmetry 
proves to also be the seed for a possible solution to this problem. 

As an example of what has just been mentioned we note from Figure (TU] 
where a fluctuation occurs just before 5000 mcs and the magnetization peaks 
at from —1. The configuration is now in the middle of the valley and happens 
to go back to its previous state. The same phenomenon occurs just after 5000 
mcs but in this instance chooses to flip to an opposite, but equally probable, 
magnetization, from -1 to 1. 

If we now were to think of the implications of this spontaneous flipping 
we come to the realization that it would cause an averaging out of the mean 
magnetization, (M). This of course has a detrimental effect on calculating the 
variance of the magnetization and thus the susceptibility. 

This can be illustrated in the Figure [T2] where the plot shows that (Af)^ 
remains zero for an extended period at low temperatures. This would cause 
the variance to peak at lower temperatures. As the lattice size increases the 
spontaneous magnetization is less likely to occur and the critical point moves 
progressively to higher temperatures, this implies that the peak for the suscepti- 
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Magnetization per spin (M/N) vs Monte Carlo steps (mcs) 
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Figure 10: This plot shows a spontaneous flip in magnetization for a 2 x 2 lattice 
at T = 1. 



bility would approach the Curie temperature from the left (lower temperatures). 
This is inconsistent with what the theory prescribes [3]. We can also conclude 
that the greater the number of mcs we use the more likely we are to introduce 
spontaneous flipping and thus averaging out of the mean magnetization which 
would move the peak of the susceptibility more to the left. 

The solution to this problem lies in a characteristic that is at the heart of 
the problem, namely that there is an equal probability for the magnetization 
to change to an opposite configuration or go back to its previous configura- 
tion. Thus if we were to use the absolute magnetization we would effectively 
be considering only one probability peak (positive one) and the complications 
of averaging out the mean magnetization would be overcome. Figure [Tl] thus 
changes to a distribution shown in Figure [T3] if we where to use the absolute 
magnetization. 

This modification doesn't come without a cost. In this instance it can be 
seen that we have averted the problems at low temperatures but end up with 
weaker statistics at high temperatures. This is apparent from the fact that 
previously we had frequent fluctuations, in Figure [TTJ between positive and 
negative magnetization at high temperatures resulting in a zero mean magne- 
tization. However we have only positive values to consider for the averaging of 
the mean magnetization producing a non zero average. This effect is reduced 
slightly since the valley is raised at high temperatures resulting in the magneti- 
zation having a high probability of being close to zero. Fortunately this nonzero 
average for the magnetization at higher temperatures is inconsequential since it 
doesn't influence the Curie temperature and appears only in the region above 
it. 

At lower temperatures the shape of the distribution changes, as indicated by 
the dotted line in Figure 1131 Thus the magnetization remains relatively stable 
in a homogeneous configuration. This is exactly the behaviour we expect and 
produces good results. 
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Figure 11: A schematic illustration of the probability density of the magneti- 
zation and how the representative spin distribution would populate a square 
latice of size L. The darker and lighter regions depict negative and positive spin 
repsectively. 



The susceptibihty that is produced in our data is thus not exactly equivalent 
to the theoretical susceptibility, x, and we will be distinguished as x'- The scal- 
ing characteristic of this susceptibility is, however, equivalent to the theoretical 
value and only varies by a constant factor above the Curie temperature. 

^ = ^ 1.1^ 22 

A comparison can be made between x' and x in Figures [14] and [15] respec- 
tively. It is clear that a marked difference in results occurs. This mistake 
becomes even more evident if you were to use x to get the critical exponent 
using finite size scaling. Only x' produces the correct finite size scaling. 

We now evaluate the plots produced using this technique discussed thus far. 
The more distinctive character of the differing plots of Figure [9] produce more 
dramatic peaks in Figure [H] of the magnetic susceptibility (x') per spin versus 
temperature as a result. This, once again, doesn't show an exact divergence 
but shows a sharp peak for the L = 16 lattice. This should be strong evidence 
eluding to a second order phase transition. 
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Figure 12: This is an illustration of the differences in the normalized values of 
(M)^; (|M|)^ and (M^) with respect to temperature. The varying lattice sizes 
considerd are 2x2 (top left); 4x4 top right; 8a;8 (bottom left) and 16a;16 (bottom 
right). 



P(M) 




+IMspl 



Figure 13: The solid line shows the revised probability density when using 
the absolute magnetization as opposed to the dotted line which represents the 
orginal propability density for magnetization. 
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Magnetic Susceptability per spin (X/N) vs Temperature (T) 
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Figure 14: This plot shows the differing results of the susceptibility for varying 
lattice sizes, L x L. 
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Figure 15: This plot shows the differing results of the susceptibility for varying 
lattice sizes, L x L. 
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2.3 Exact Results 

As pointed out in the beginning of this work there are significant chaUenges to 
the calculation of exact solutions but we do need to evaluate the dependability 
of the numerical process used in this discussion. For this purpose we restrict 
the comparison between simulation and exact results to only a 2 x 2 This only 
serves as a good first order approximation but the corroboartion should imporve 
as the lattice size is increased. 

The 2** different configurations that can be listed can be reduced by sym- 
metric considerations to only four configurations. Configuration A is the fully 





Figure 16: The four significant lattice configurations for a 2 x 2 lattice. 

aligned spin configuration and configuration B has one spin in the opposite di- 
rection to the other three while C and D have two spins up and two spins down. 
To generate the full sixteen different configurations we need only consider geo- 
metric variations of these four configurations. 

Using equations (118p and p7p we can calculate the energy and magnetization 
for each of the four configurations. Taking into account the degeneracy of each 
of these configurations allows us to generate the exact equations for calculating 
the relevant observables. These results are listed in Table [TJ 



Configuration 


Degeneracy 


Energy(£') 


Magnetization(M" 


A 


2 


-8 


+4,-4 


B 


8 





+2,-2 


C 


4 








D 


2 


8 






Table 1: Energy and Magnetization for respective configurations illustrated in 
Figure [TH 



Using equations ^ and © in conjunction with the results of Table [T] pro- 
duces the following equations 
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Z = 2 e^^-^ + 12 + 2 e-*'' (23) 

(i;) = -^[2(8)e«'^ + 2(-8)e-8/3] (24) 

= \ [ 2(64) e^'^ + 2(64) e"**'' ] (25) 

(|M|> = i[2(4)e«'3 + 8(2)] (26) 

(Af2) = 1[ 2(16) e«^ + 8(4)]. (27) 

Zj 



Applying equations ([21]); ([251); ([Ml) and ([27l) to the formulas given in (fTg]) 
and p2|) we can obtain the heat capacity and susceptibility respectively. The 
comparison of these results are shown in Figure [T71 The results achieved by the 
Monte Carlo method match the exact calculation exceptionally well and only 
deviate very slightly from the prescribed heat capacity at high temperatures. 



Specific Heat Capacity per spin (C/N) vs Temperature (T) 
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Susceptibility per spin (X'/N) vs Temperature (T) 
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Figure 17: This plot shows a favourable comparison between the exact calcula- 
tions done for the Heat Capacity per spin (C /A^ ) and the Susceptibility per spin 
(x'/A^) with the numerically generated solutions of the Monte Carlo simulation. 



2.4 Finite size scaling 

One of the limitations that the Ising Model confronts us with is the finite size of 
our lattice. This results in a problem of recognizing the specific point at which 
the phase transition occurs. This should be at a theoretical point of divergence 
but we are limited by the size of the lattice under consideration and thus dont see 
this divergence. This effect is minimized by using periodic boundary conditions 
but would only be resolved if we where to consider an infinitely sized lattice 
as with the associated theoretical values for the phase transition. It is thus 
necessary to use a construct that will allow us to extrapolate the respective 
theoretical value given the limited resource of a finite sized lattice. The aptly 
named procedure of finite size scaling is used to do just this. 

It becomes useful to define a critical exponent to better understand the 
nature of the divergence near the critical temperature. The critical exponent, 
A, is given by A = limt^o [fj^ ^ or more commonly written as F(t) ~ \t\^ 
where t — (T ~ Tc). This exponent is important in the sense that it offers a 
more universal characteristic for differing data collected. This attribute will be 
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taken advantage of and illustrated by showing that in a reduced unit plot the 
data collapses to a common curve with the same critical exponent. 
The critical exponents relevant to the Ising model are as follows: 

M{T) - (Tc - Tf 
C IT-T^r" 



If we examine the relationship between the lattice size, L, with respect to the 
temperature relationship, |T — Tc\, we discover that \T ^ Tc\ « 1 asL^oo. 
Thus a critical exponent is also applicable for the lattice size. This produces 
L ~ \Tc{L — oo) — Tc{L)\~^ which can in turn be used to reduce the above 
exponents for the Ising model to a more appropriate form, in terms of lattice 
size. 



£^{T)^\T-T,\-'' (28) 



M(T) - 


' (T, - Tf L-P'" 


(29) 






(30) 


X ' 




(31) 



This proportionality is illustrated in a schematic form in Figure [TH] of the re- 
lationship between the plot of a certain observable and the respective critical 
exponent. In this instance susceptibility was used as an example. 




L 



Figure 18: Theoretical illustration of Finite Size Scaling relationship for suscep- 
tibility X- 

We can now attempt to use equations ((25)) - (PT|) and determine their ap- 
propriate exponents. This is simply done by taking the peak values for the 
collected data of the observables and plotting a In-ln graph that should yield a 
straight line with the gradient being equal to the respective critical exponents. 
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This procedure is made easier since u is equal to 1 for a two dimensional lattice. 
An additional calculation for the observables of a L = 32 lattice were done, 
since this would increase the number of data points. This would also offset the 
poor statistics associated with the L — 2 lattice and ultimately allow for a more 
accurate result. 



Ln of Susceptiblity (In X) vs Ln of Lattice size (In L) 
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Figure 19: Plots obtained for the calculation of the critical exponents. 

In Figure [19] the calculation for the critical exponent is shown. The answers 
are listed in Table [H The graph for the heat capacity isn't a straight line and 
shows a curvature. The reason is that a is zero in the two dimensional Ising 
model and should rather be interpreted as C ~ Cq In L. 



Quantity 


Exponent 


Finite Size 
Scaling (2d) 


Theoretical (2d) 


Magnetization 


P 


0.128 ± 0.008 


0.125 


Susceptibility 


7 


1.76 ± 0.01 


1.75 


Heat Capacity 


Co 


0.518 ± 0.02 


0.500 



Table 2: Displays the calculated and theoretical critical exponents. 



We can in a naive sense attempt to determine the Curie temperature by 
implementing the critical relation for the lattice size. This is listed in Table [31 
It doesn't prove to be very accurate and we thus have to implement a different 
notion in order to ascertain the critical temperature accurately. 

The transition point can be determined by using the cuniulant 



(32) 
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2(T-Tc) vs (X/N)*2 
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Figure 20: This graph shows a reduced unit plot for x' with a critical exponent 
of 7 = 1.75. 



Lattice Size (L) 


Estimated Tc 


Estimated 


Theoretical 






Tc{L oo) 


TciL = oo) 


2 


3.0 


2.50 


2.269 


4 


2.8 


2.55 


2.269 


8 


2.5 


2.43 


2.269 


16 


2.4 


2.34 


2.269 


32 


2.3 


2.27 


2.269 



Table 3: Listed critical temperatures calculated from the critical lattice size 
relation. 



This calculation has to use double precision in order to retain a high level of 
accuracy in calculating the critical point. To determine the critical point we 
choose pairs of linear system sizes (L, L'). The critical point is fixed at /7l = Ul' ■ 
Thus taking the ratio of different cumulants for different sized lattices, Ul/Ul', 
we will get an intersection at a particular temperature. This is the desired 
critical temperature. This procedure is not as straightforward as it may seem 
and requires the cumulants to be collected very near to the transition point. 
Thus an iterative process may need to be employed in order to narrow down 
the region of where the critical temperature is located. 

This analysis is done until an unique intersection is found for all cumulants. 
This method is illustrated in Figure HH The L = 2 lattice isn't shown since it 
doesn't exhibit the level of accuracy that we desire. From the fitted graph of 
the cumulants it can be seen that intersection is common to all the cumulants 
except for the U4. This points towards a poor value for this particular statistic 
and is thus not used, along with U2, the ratios of the cumulants to determine 
the Curie temperatures. The final analysis produces a result of a Tc = 2.266 
which agrees favourably with the theoretical value of Tc ~ 2.269. 
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Figure 21: Plots of the cumulants {Ul)- 



2.5 Conclusion 

In this review we pointed out why an exact exposition of an Ising model for 
a ferromagnet is not easily achieved. An argument was proposed and moti- 
vated for making the computational problem far more tractable by considering 
a stochastic process in combination with the Metroplis-Hastings sampling algo- 
rithm. The numcric;al results produced by the Monte Carlo simulation compare 
favourably with the theoretical results and are a viable and efficient alternative 
to an exact calculation. The requirements for producing accurate results are 
to consider large lattice sizes and a large number of Monte Carlo steps. The 
accuracy is very compelling even for small lattice sizes. 

It is important to note and make provision for the potential of spontaneous 
magnetization to occur in a finite sized lattice. This can have serious conse- 
quences on the accuracy of the magnetization and the susceptibility which in 
turn will lead to incorrect results for the finite size scaling of these observables. 
The occurrence and severity of spontaneous magnetization is directly propor- 
tional to the number of Monte Carlo steps used and inversely proportional to the 
lattice size considered. A practical means to overcome this complication is to use 
the absolute magnetization in the variance of the susceptibility instead of just 
the magnetization. This is an effective solution that produces good statistics 
only deviating slightly at high temperatures from the theoretical values. 

In conclusion a finite size scaling analysis was undertaken to determine the 
critical exponents for the observables. These where in good agreement with 
the theoretical exponents. The critical temperature was also calculated using a 
ratio of cumulants with differing lattice sizes and generated results which were 
in good agreement with the theoretical values. 
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A Monte Carlo source code 



Monte Carlo Simulation of a 2D Ising Model Page 1/2 



#lncliide <iostream . h> 

#±nclude <m.ath , h> 

#lnclude <stdlib.h> 

#lnclude <f stream . h> 



from "Numerical Recipes in C" (Rani , c) 



//file to output data into 

of stream DATA ( "DATA.l.dat" , ios : :out); 



//structure for a 2d lattice with coordinates x and y 

struct lat_type 

{ 

int x; 
Int y; 

}; 



-1; 



floBt T=5 . 0; 

const float minT=0.5; 

float change=0 . 1; 

int lat [size + 1] [size + 1] ; 

long unsigned int mcs=10000; 

int transient=1000; 

double norr[i= (1 , 0/float (mcs*n} ) ; 

long int seed=436675; 



//lattice size 
/ /array size for lattice 
//number of spin points on lattice 
//starting point for temperature 
//minimum temperature 

//size of steps for temperature loop 
//2d lattice for spins 
//number of Monte Carlo steps 
//number of transient steps 
//normalization for averaging 
//seed for random number generator 



//function for random initialization of lattice 
initialize (int lat[size+l] [size+1]) 

{ 

foE(int y=size;y>=l; y — ) 

:=size;x++) 



for (ifit 

{ 



if (rani (sseed) >=0 . 5) 
lat[x] [y]=1; 

else 

lat[x] [y]=-l; 



'output of lattice configuration to the screen 
tput(ir,t lat [size+1] [size+1] ) 



forfirtt x=l;x<=size;x++) 
{ 

if (lat[x] [y]<0) 

cout<<" — " ; 
else 

cout<<" + " ; 

} 

cout<<endl; 



'/function for choosing random position on lattice 
hoose_random_pos_lat ( lat_type spos ) 

pos . x= (i/)t) ceil (rani (sseed) * (size) ) ; 
pos .y= (int) ceil (rani (sseed) * (size) ) ; 
if (pos , x>size I |pos.y>size) 
{ 

cout<<"error in array size allocation for random position on lattice!"; 

exit; 

1 



//function for calculating energy at a particular position on lattice 

int energy_pos (lat_type Spos) 

{ 

//periodic boundary conditions 
int up, down , lef t , right , e ; 
if (pos . y==size) 
up = l; 

else 

up=pos.y+l; 
if (pos . y==l) 

else 

down=pos . y-l ; 
if (pos . x==l ) 

left=size; 

else 

left=pos . x-1; 
if (pes . x==size) 
right =1; 

else 

right =pos -x+l ; 

//energy for specific position 
e=-l*lat [pos -x] [pos.y] 

* (lat [left] [pos .y]+lat [right] [pos .y] +lat [pos.x] [up] +lat [pos.x] [down] ) ; 
return e ; 



type pos 
s (pos) ; 



else if (rani (S'-;eed) <exp (-de/T) } 
else 

return false; 



//change 1 

//flip due 
//flip due 
//no flip 



I energy for specific spin 
to Jopver energy 
to heat bath 



//flip spin at give 
f lip (lat_type pos) 



lat [pos .x] [pos . y] =-lat [pos . x] [pos . y] ; 
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//function for disregarding transient results 

trans lent _results ( ) 

( 

lat_type pos; 
int de=0; 

for (int a=l; a<=transient; a++> 
{ 

for(int b=l;b<=n;b++) 
{ 

choose_random__pos_lat (pos) ; 
if (test_f lip (pos, de) ) 



//function for calculating total magnetization of lattice 
int total_magnetization () 

{ 

i;;r n 0; 

for(^;.; V size;y> = l;Y — ) 

for(int x=l;x<=size;x++) 
( 

m+=lat[x] [y]; 

} 



//function for calculating total energy of lattice 

int total_energy ( ) 

( 

lat_type pos; 
int e=0; 

for (int y=slze; y>=l; y — ) 



return e; 



//declaring- variables to be used in calculating the observables 

double E=0, Esq=0, Esq_avg=0 , E_avg=0 , etot=0,etotsq=0; 

double M=0 , Msq=0 , Msq_avg=0 , M_avg=0 , rrLtot=0 , mtotsq=0 ; 

double Mabs=0, Mabs_avg=0,Mq_avg=0,mabstot=0,mqtot=0; 

int de=0; 

lat_type pos; 



/ /Tempera t ure 1 oop 

for (; T>=minT; T=T-change) 



s () ; 



ns values 



//observables adopt equilibrated lattice conflgurati 
M=total_magnetization {) ; 
Mabs=abs (total_magnetization ( ) ) ; 
E=total_energy () ; 

//initialize summation variables at each temperature step 

etot=0; 

etotsq=0; 

mtot=0; 

mtotsq=0; 

mabstot=0; 

mqtot=0; 



choose_random_pos_lat (pos) ; 

if (test_f lip (pos, de) ) 

{ 

flip (pos) ; 

//adjust observables 
E+=2*de; 

M+=2*lat [pos.x] [pos -y] ; 
Mabs+=abs (lat [pos ,x] [pos . y] ) ; 

} 

) 

//keep summation of observables 

etot+=E/2 . 0; //so as not to count the energy for each spin twice 

etotsq+=E/2 . 0*E/2 .0; 

mtot+=M; 

mtOtsq+=M*M; 

mqtOt+=M*M*M*M; 

mabstot+= (sqrt (M*M) } ; 

} 

/ /average observables 
E_avg = etot *norm; 
Esq_avg=etotsq*norm; 
M_avg=mtot *norm; 
Msq_avg=mtotsq* norm; 
Mabs_avg=mabstot *norm; 



Mq_ 



mqtot* 



//output data to file 

DATA<<T<< 

"\t"<<M_avg<<"\t"<<Mabs_avg<<"\t"<<Msq_avg<< 

"\t"<< (Msq_avg- (M_avg*M_avg*n) ) / (T) << 

"\t"<< (Msq_avg- (Mabs_avg*Mabs_avg*n) ) / (T) < 

"\l"<<E_avg<<"\t"<<Esq_avg<< 

"\t"<< (Esq_avg- (E_avg*E_avg*n) ) / (T*T) << 

"\t"«l- ( (Mq_avg) / (3*Msq_avg) ) «endl; 



//t empera t ure 

//<M>; < I Ml >; <M''2> per spin 
//susceptibility per spin (X) 
'•//susceptibility per spin (X' ) 
//<E>; <E"2> per spin 
//heat capacity (C) per spin 
//cumulant (U_L) 
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